/** \file    axis.h
  * \brief   Axis = directionless vector: with unit length or not.
  * \date    2002
  * \author  Pavel Agapov, Fedor Chelnokov
  * \version $Header: /home/cscom1r4/cvsroot/cs/include/typelib/axis.h,v 1.1 2004/09/21 18:21:30 cscom1r4 Exp $
  */

#ifndef _axis_h_
#define _axis_h_

#include "typelib/types.h"
#include <limits>

namespace types {

template <class T>
class axis : public point<T> { 
public:
	explicit axis(T x = 0, T y = 0)
		: point<T>(x, y) { }
	explicit axis(point<T> p)
		: point<T>(p) { }
	axis & normalize() {
		double a = _hypot(x, y);
		if (a > std::numeric_limits<T>::epsilon())
			(* this) /= (T)a;
		else
			x = y = 0;
		return * this;
	}
	axis & not_negative_x() { 
		if (x < 0) 
			inverse(*this); 
		return *this;
	}

	axis<T> & operator += (const axis<T> & a) {
		point<T> p(a);
		if (scalar_product(*this, a) < 0)
			inverse(p);
		static_cast<point<T> &>(*this) += p;
		return *this;
	}
	axis<T> & operator -= (const axis<T> & a) {
		point<T> p(a);
		if (scalar_product(*this, a) < 0)
			inverse(p);
		static_cast<point<T> &>(*this) -= p;
		return *this;
	}
	friend axis<T> operator + (const axis<T> & x, const axis<T> & y)
		{ axis<T> t(x); return t += y; }
	friend axis<T> operator - (const axis<T> & x, const axis<T> & y)
		{ axis<T> t(x); return t -= y; }
};

} //namespace types

template <class T>
inline void inverse(axis<T> & a) 
	{ inverse(static_cast<point<T> &>(a)); }

template <class T>
inline double scalar_product(const axis<T> & a0, const axis<T> & a1) {
	const point<T> & p0(a0), & p1(a1);
	return scalar_product(p0, p1);
}

//!
//! general rule: normalization should be called by user explicitly
//!

template <class T, class F>
inline axis<T> linear(const axis<T> & a0, const axis<T> & a1, F a) {
	point<T> p0(a0), p1(a1);
	if (scalar_product(p0, p1) < 0)
		inverse(p1);
	return axis<T>(linear(p0, p1, a));
}

template <class T, size_t N>
inline void adjust_dir(point<T> *p) {
	double max = 0;
	int i, j, index;
	for (i = 0; i < N; ++i) {
		double v = 0;
		for (j = 0; j < N; ++j)
			v += fabs(scalar_product(p[i], p[j]));
		if (v > max) {
			index = i;
			max = v;
		}
	}
	if (max == 0)
		return; // axis<T>();
	for (j = 0; j < N; ++j) {
		if (index != j && scalar_product(p[index], p[j]) < 0)
			inverse(p[j]);
	}
}

template <class T, class F> 
inline axis<T> linear(const axis<T> & x0, const axis<T> & x1, const axis<T> & x2, F a1, F a2) {
	point<T> p[3] = {x0, x1, x2};
	adjust_dir<T,3>(p);
	return axis<T>(linear(p[0], p[1], p[2], a1, a2));
}

template <class T, class F>
inline axis<T> bilinear(const axis<T> & a00, const axis<T> & a01, const axis<T> & a10, const axis<T> & a11, F a0, F a1) {
	point<T> p[4] = {a00, a01, a10, a11};
	adjust_dir<T,4>(p);
	return axis<T>(bilinear(p[0], p[1], p[2], p[3], a0, a1));
}

#endif //_axis_h_
